Connecting systems with short and long ranged interactions: local molecular field 

theory for ionic fiuids 
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Structural and thermodynamic properties of ionic fiuids are related to those of a simpler "mimic" 
system with short ranged intermolecular interactions in a spatially varying effective field by use of 
Local Molecular Field (LMF) theory, already successfully applied to nonuniform simple fluids. By 
consistently using the LMF approximation to describe only the slowly varying part of the Coulomb 
interaction, which we view as arising from a rigid Gaussian charge distribution with an appropriately 
chosen width a, exceptionally accurate results can be found. In this paper we study a uniform system 
of charged hard spheres in a uniform neutralizing background, where these ideas can be presented 
in their simplest form. At low densities the LMF theory reduces to a generalized version of the 
Poisson-Boltzmann approximation, but the predicted structure factor satisfies the exact Stillinger- 
Lovett moment conditions, and with optimal choice of a the lowest order approximation remains 
accurate for much stronger couplings. At high density and strong couplings the pair correlation 
function in the uniform mimic system with short ranged interactions is very similar to that of the 
full ionic system. A simple analytic formula can then describe the difference in internal energy 
between the ionic system and the associated mimic system. 



I. INTRODUCTION 

In this paper we describe a new theory for the structure 
and thermodynamics of ionic fluids based on a general- 
ization of the local molecular field (LMF) theory we have 
successfully applied to nonuniform simple fluidsi ^i^i'^i'^i^i^ 
A basic step in the LMF theory is the replacement of 
longer ranged and slowly varying parts of the intermolec- 
ular interactions by an appropriately chosen effective sin- 
gle particle potential. The structure and thermodynam- 
ics of the resulting reference system with shorter ranged 
intermolecular interactions in the presence of the effec- 
tive field is then related to that of the original system of 
interest. 

This strategy seems particularly appropriate for ionic 
systems since at long distances the Coulomb interaction 
is weak and very slowly varying, and systems with short 
ranged interactions are significantly easier to treat by 
theory or simulations. However Coulomb interactions 
can be very strong and rapidly varying at short distances. 
A key question we address is how to divide the Coulomb 
interaction into "short" and "long" ranged parts so that 
the LMF theory can give accurate results. Its answer 
allows for the first controlled use of the LMF approxima- 
tion, and we find exceptionally accurate results, even bet- 
ter than those found earlier for fiuids with short ranged 
interactions. 

Although the most physically interesting applications 
of these ideas are probably to nonuniform mixtures of size 
and charge asymmetric ions, in this initial discussion we 
consider a uniform one component charged hard sphere 
system (OCCHS) where almost all the ideas in the LMF 
theory can be seen in their simplest formii The OCCHS 
is made up of (say positively) charged "ions" comprised 



of hard spheres with a diameter d > with embedded 
positive point charges in the presence of a uniform neu- 
tralizing background. The only nontrivial correlations 
are between the positive ions and for most purposes we 
can think of this as a one component system with very 
long ranged repulsive interactions. A special case is the 
one component plasma (OCP) where there is no hard 
core (d = 0). Nothing in the theory makes essential use 
of the simplifications in the OCCHS. Generalizations to 
nonuniform and asymmetric models are straightforward 
in most cases, and equally good results have been found. 



II. LOCAL MOLECULAR FIELD EQUATION 
A. Nonuniform systems 

We first discuss the qualitative ideas leading to the 
LMF equation. These will be further developed and 
made more precise in our discussion of the OCCHS. The 
simplest application of LMF theory relates the structure 
and thermodynamics of a nonuniform system of interest 
with a spherically symmetric pair potential w{r) in an ex- 
ternal field whose value at any point ri is (piri) to those 
of a reference system with a shorter ranged pair interac- 
tion uo{r) in a renormalized effective field (/</j(ri). (jjji is 
supposed to be chosen to take account of the averaged 
effects of the perturbation potential ui{r), where 



w{r) — Uo{r) + Ui{r). 



(1) 



This separation of the intermolecular interaction w into 
two parts can be done in an infinite number of ways, and 
any choice of ui defines a possible associated reference 
system. However the averaging procedure leading to the 



2 



Potential Separation 




r 



Figure 1: Tlie separation of tlie 1/r potential into a short 
ranged piece uo{r) = eTic{r/a)/r and a long ranged piece 
ui(r) = eri{r/a)/r. A bigger a corresponds to a longer ranged 
mimic system uo(r), and a more slowly varying ui{r). Here 
two relevant a values are shown for comparison. 

simple LMF theory can be expected to give very accurate 
results only for certain properly chosen slowly varying ui. 

Figure 1 gives examples of separations of the repulsive 
Coulomb potential we will use in this paper, parameter- 
ized by an important length scale a. As explained in 
detail in Section llVI when a is chosen larger than some 
state-dependent minimum size CTmin, the Coulomb per- 
turbation ui is sufficiently slowly varying that the LMF 
theory can give very accurate results. This is the crucial 
step in developing a simple and accurate theory for ionic 
systems. 

We will refer to the resulting special reference systems 
with properly chosen a as "mimic systems." In realistic 
models of ionic solutions there are always strong short 
ranged repulsive core interactions that must be dealt with 
in any quantitative theory or simulation. The mimic sys- 
tem simply treats the short ranged rapidly varying part 
of the Coulomb potential as an additional core-like con- 
tribution that generates a modified core interaction. As 
we will see, many properties of the full long ranged sys- 
tem can be very accurately described using those of the 
short ranged mimic system. 

For simple fluids with short ranged interactions the 
LMF approach has proved most useful when w can be di- 
vided into a slowly varying perturbation ui{r) describing 
the relatively weak and longer ranged attractive interac- 
tions and a short ranged rapidly varying core potential 
Mo, which accounts for the local excluded volume correla- 
tions of the particles.? A separation with these qualitative 
features suffices to motivate the derivation of the basic 
LMF equation that follows. 

For any given (j){r) an associated 0fl(r) could always 
found in principle so that the nonuniform singlet density 
po(r; in the reference system (denoted by the sub- 
script 0) equals that in the full system p(r; [0]). Of course 
the latter is not known in advance and its determination 



is one of the main goals of the theory. However, if a ui (r) 
can be chosen to be slowly varying over the range of ex- 
cluded volume correlations induced by the short ranged 
potential uo{r), then we can make some physically moti- 
vated approximations to derive a self-consistent equation 
to determine the associated (/jr. 

At high densities we expect that short ranged corre- 
lations in both systems are controlled by packing effects 
from the identical repulsive cores, and it seems plausible 
that can be chosen so that both the singlet densities 
and the conditional singlet densities in the reference and 
full systems resemble each other. That is, when 

po{r;[cl>R]) = p{r;[cl>]), (2) 

we also expect that 

/7o(ri|r2; [(/>/?]) ~ p(ri|r2; [(/)]) (3) 

holds to a good approximation. Here po(ri|r2; [(t^R]) is 
the (conditional) density at ri given that a particle is 
fixed at r2, directly related to the nonuniform pair cor- 
relation function. With this assumption we can derive 
an equation for (j)R that also turns out to give exact re- 
sults at very low densities, where pair correlations are 
not important. 

As discussed previously,'* by subtracting the balance 
of forces as described by the exact Yvon-Born-Green 
hierarchjii for the full and reference systems we find a 
relation between the associated forces 

-Vi[(?!)i?(ri)-(/)(ri)] - j dr2Po{Y2\Yi\ ['/'/?])ViUi(ri2). 

(4) 

Moreover, if ui(ri2) is very slowly varying over the 
range of short ranged pair correlations, then ViUi(ri2) 
essentially vanishes in the range of integration where 
po(r2|ri; [<I>r]) differs significantly from yOo(r2; [(I>r]) in eq 
0] Then we can replace the former by the latter, and take 
the gradient outside the integral and integrate eq^J 

In making this replacement we have ignored correla- 
tions between the particles at ri and r2, and in most 
contexts this would represent a crude and generally inac- 
curate approximation. However for slowly varying ui we 
see that this particular use of the (mean field) approxima- 
tion can be very accurate, even at high density. Choosing 
the constant of integration so that the bulk densities in 
zero field satisfy = Pq we arrive at the simple local 
molecular field (LMF) equationii^ for the effective field 
4>R- 

(t)R{ri) = (t)iri) + J dr2[po{r2; [(t)R]) ~ Po]ui{r 12), (5) 

which is the starting point for our work on nonuniform 
fluids with both short and long ranged interactions. 

To solve this self consistent equation we need to de- 
termine the nonuniform density po(r2; [4'r]) in the pres- 
ence of the effective field cfiR. The LMF approach does 
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not specify or require a particular way to do this. How- 
ever since the intermolecular interactions and the effec- 
tive field have shorter ranges in the reference system, 
both theory and simulations of the nonuniform structure 
are usually easier to carry out than in the full system. 
Equation 121 then provides the fundamental link between 
structure in the nonuniform reference and full systems, 
and from this thermodynamic functions can be deter- 
mined. 

The name "local molecular field" is suggested by the 
direct analogy to the spatially varying effective field in- 
troduced in the usual mean or molecular field theory for 
a nonuniform Ising modelAiSi However, while the lat- 
ter theory is usually viewed as a crude approximation, 
the derivation sketched above suggests that if a proper 
choice of a slowly-varying ui is made, then accurate re- 
sults should be found from a self-consistent solution of 
the LMF equation in many cases, provided that an ac- 
curate treatment of the density response po(v; [(j)R]) to a 
given is used. In this paper we use the exact eqlHlbe- 
low at low densities and results of computer simulations 
at higher densities, so whatever errors remain arise only 
from the LMF approximation itself. 

B. Uniform systems 

LMF theory is equally useful for uniform fluidsiS*ii In 
particular eq[5lis consistent with the physical idea that 
in a dense uniform fluid with (j) — the forces associated 
with the slowly varying ui from oppositely situated parti- 
cles essentially cancel in most relevant configurationsi^ii^ 
Moreover, any residual effects are strongly damped by the 
small compressibility at high density, so we expect that 
the radial distribution functions will satisfy 

5o(r) ~ gir), (6) 

as predicted by eqs |21 and [S] for (f) = (pji — 0. For exam- 
ple, eq El holds to a rather good approximation at high 
density in the uniform Lennard- Jones (LJ) fluid provided 
that the WCA separation2ii2ii4 with its relatively slowly 
varying ui is used, showing the consistency of the physi- 
cal picture. 

This perfect cancellation argument can at best give 
reasonable results only for uniform fluids at high den- 
sity. However, LMF theory can be applied to a general 
nonuniform fluid, and by taking such a perspective and 
using only the basic eq [3 we can significantly improve 
on the predictions of eq^for pair correlations in uniform 
fluids.^'" 

Corrections to eq El can be found by considering the 
particular external field arising from a fluid particle fixed 
at the origin, </)(r) = w{r). The induced density now gives 

p(r;H)=p^g(r). (7) 

This exact equationi^ relating the nonuniform singlet 
density induced by a fixed particle to the radial distri- 
bution function g(r) in the uniform fluid plays a key role 



in the theory below. There are now net unbalanced forces 
arising from the fixed particle and eq[31predicts a nonzero 
(pR, which can be used in eqOto give a more accurate 
approximation for p(r; [w]) and hence g{r). 

For the uniform LJ fluid this approach accurately de- 
termines the small corrections to eqjSlat high densityi^iii 
Moreover at very low densities where eqjHlwould be very 
inaccurate, eqjSlgives (/li? (r) = 0(r) = w{r) and we obtain 
the exact low density limilii for g{r) = exp[— /3w(r)] by 
using eq[21and the exact low density limit of the reference 
system in the field (/jr: 

Po(r;[0fl])=p^e-'='*-('^). (8) 

Here f3 ^ (kBT)-^. 

Similar accurate results for the nonuniform LJ fluid 
have been found for more general external fields repre- 
senting hard core solutes of various sizes, for the liquid- 
vapor interface, and for drying transitionsi^ Thus LMF 
theory has provided a qualitatively and often quantita- 
tively accurate description of structure, thermodynam- 
ics, and phase transitions in fluids with short ranged 
interactions >2i 



III. ONE COMPONENT CHARGED HARD 
SPHERES 

In this paper we focus on the uniform OCCHS, where 
there are N positive ions in a volume V with a uniform 
neutralizing background that also penetrates the ions. 
The pair potential ■w{r) for the ions in the OCCHS is 
usually written as 

w{r) = Wdir) + Wq{r), (9) 

the sum of a hard sphere potential 

and the pair potential Wq{r) arising from point charges 
of magnitude q, where 

Wgir) = ^. (11) 

The separation in eqjSlis a special case of eq^and more 
general separations of w will prove useful in the LMF 
theory developed below. In eq^the solvent is crudely 
represented by a uniform static dielectric constant e. In 
the hmit d = 0, the OCCHS reduces to the OCP. 

It is convenient to introduce a characteristic length de- 
scribing the typical distance between neighboring parti- 
cles. A standard choice is the ion sphere radius a chosen 
so that 



The nearest neighbor spacing is about 1.6a when the ions 
are arranged in a simple cubic lattice. 
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Thermodynamic properties in the OCCHS can then 
be characterized using two dimensionless parameters, the 
ionic strength 

T.^A (13) 
ea 

which compares the bare Coulomb interaction energy be- 
tween two ions separated by the characteristic distance a 
to fcsT, and the hard sphere packing fraction 

ri = (14) 

Note that d/a = 2rj^/^. In the OCP d = 77 = and 
thermodynamic properties depend only on the single di- 
mensionless parameter F. 

Pair correlations between the ions in the uniform fluid 
are most conveniently described in terms of the density 
change induced by fixing a particle at the origin: 

Ap(r; H) = p{r; [w]) - = p^/i(r), (15) 

where = N/V and h{r) = g{r) — 1 is the pair correla- 
tion function in the uniform fluid. 

The unique consequences of the long ranged interac- 
tion in the OCCHS are most easily seen by taking the 
Fourier transform of eq 1151 and defining the dimension- 
less structure factor 

S{k) = l + /:ip{k-[w]) = l^ p^h{k), (16) 

where the caret denotes a Fourier transform. As argued 
generally in the seminal work of Stillinger and Lovett^^ 
(SL), there should be complete screening at long wave- 
lengths of any induced charge distribution in a conduct- 
ing ionic fluid. This constrains the behavior at small 
wavevectors of the charge-charge correlation function. 
For the OCCHS the only nontrivial correlations are be- 
tween the positive ions and the results of SL reduce to 
the requirement that S{k) has the universal form 

s{k) = {) + e/ki + o{k''), (17) 

independent of any details of the short ranged core po- 
tential Wd or any other short ranged interactions that 
might exist. Here k^ is the Debye wavevector, deflned 

by 

fc|) EE 47r/3gV^/e = Sr/a^. (18) 

The exact vanishing of S'(fc) at fc = arises from elec- 
trical neutrality (the "zeroth" moment condition) and 
the fixed coefficient of the quadratic term is an example 
of the famous SL second moment conditioni^ii^ This be- 
havior is distinctly different than that found in any fiuid 
with short ranged interactions, where So{k) at fc = is 
finite, proportional to the compressibility, and the coef- 
ficient of k'^ depends on the details of the intcrmolccular 
interactions. 

Thermodynamic properties can be found by integra- 
tion of the correlation functions. In particular, account- 
ing for the background by taking the appropriate limits 



of the standard result for a two component system^iS, the 
excess internal energy (over the ideal gas) of the uniform 
OCCHS can be exactly written as 

f3E-- PpB f q\ 



IV. LMF THEORY FOR OCCHS 
A. Gaussian charge distribution 

We want to use the general LMF equation to describe 
the ion correlation functions in the uniform OCCHS. It 
is clear from the derivation in section ^] that a proper 
separation of the interaction potential w ~ uq + ui is 
required for this self-consistent approach to be accurate. 
At first glance it may seem natural to use the separation 
on the right side of eq IHl where uq is taken to be the 
hard core potential Wd and ui is the full point charge 
interaction Wg in eg 1111 However at short distances out- 
side the hard core the Coulomb potential can be strong 
and rapidly varying and such interactions cannot be ac- 
curately treated by the averaging used in eq[5l 

Indeed in the OCP with no hard core there are arbi- 
trarily large and rapidly varying interactions as r —^ 0. 
This limit makes it clear that we should try to separate 
the point charge Coulomb pair interaction Wq{r) itself 
into a slowly varying part Wqi{r), which we will take as a 
particularly appropriate ui{r) to use in the LMF theory, 
and then combine the remainder WqQ{r) = Wq(r) — ui{r) 
with Wd{r) (and more generally with any other short 
ranged core interactions that exist) to give the associ- 
ated uo{r). 

This strategy differs from that used in many density 
functional and integral equation methods, where one first 
chooses a mathematically convenient or especially simple 
reference potential uo{r) and then treats the remainder 
w{r) — uo{r) as a perturbation, taking advantage of the 
particular form of the reference system in whatever ap- 
proximate theories are used. However, for Coulomb sys- 
tems at least, the existing theories often have large and 
uncontrolled errors with the usual choices of uq. We be- 
lieve our approach of choosing a slowly varying ui for use 
in the LMF theory offers many conceptual and compu- 
tational advantages, and it connects directly to similar 
physically motivated work on fluids with short ranged 
interactions. 

To that end we interpret the 1/r term in Wq as the elec- 
trostatic potential of a unit point charge in vacuum. The 
same slowly varying asymptotic behavior would come 
from any other normalized charge distribution and the 
"smearing" of the point charge would produce a less 
rapidly varying potential at small r. This suggests us- 
ing a properly chosen charge distribution to determine 
Wqi{r) or Mi(r). 

Consider in particular as in the Ewald sum method^iiS 
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a normalized unit Gaussian charge distribution 

Pair) = exp[-(r/a)2], (20) 

where a indicates the length scale of the smearing. The 
particular advantages of this choice will soon become ap- 
parent. Our use of a Gaussian charge distribution to 
determine a slowly varying part of the Coulomb pair in- 
teraction is simpler than in the Ewald sum method, which 
considers periodic images of ion configurations with em- 
bedded screening (negative) and compensating (positive) 
Gaussian charge distributions. -'■^ 

The electrostatic potential Va- (r) arising from eal20lsat- 
isfies Poisson's equation V^W(j(r) = ~~4:TTPa-{r), which is 
easily solved by Fourier transform to give 



or m r-space, 



47r 



Va{r) = - erf(r/(T), 



(21) 



(22) 



where erf is the usual error function. The point charge 
model corresponds to the limit (7 = 0. 
Thus we can write 

i = ierfc(r/(7) + - erf(r/(7), (23) 
r r r 

and use this identity to extract from the dimensionless 
Coulomb pair interaction (3wq{r) a cr-dependent pertur- 
bation piece (3ui{r): 

j3ui{r) = —v,{r) = ^ erf(r/a). (24) 
e er 

This perturbation remains finite as r — > 0, with /3ui(0) = 
27r-i/2/3g2/(g^)^ 

As illustrated in Figure ^ with appropriate choices of 
a we can produce a very slowly varying ui{r), which from 
eg 1211 also decays very rapidly in fc-space. As argued in 
sectional these are the qualitative features that would be 
most appropriate for a perturbation ui to give accurate 
results from LMF theory. The choice of a in the Gaus- 
sian charge distribution permits a controlled use of the 
LMF approximation, and as shown below, with proper 
choices of a the LMF theory can give exceptionally ac- 
curate results. 

Any particular choice of a in eql^then fixes the as- 
sociated reference system interactiorkSS as 



/3uo(r) = I3ud{r) + 



rfc(r/CT). 



(25) 



The Coulomb part w^o {r) of the reference interaction de- 
cays very rapidly for r > a. For large r we have 



(3wqo{r) = — -erfc(r/CT) 



■exp[-(r/a)2]. (26) 



We call the special reference systems that result from 
optimal choices of a as discussed in section IIV Dl below 
mimic systems, since at high density the local structure 
in the uniform mimic system as exhibited in gi){r) very 
accurately approximates the g{r) of the full system as in 
eqEl This again illustrates the consistency and accuracy 
of the LMF approach when an appropriate mimic system 
is used. 

Equation 121 can also be interpreted physically as the 
Coulomb energy arising from two ions each with a rigid 
Gaussian charge distribution, eq 1201 with a width a — 
a l\pl. More generally, in ionic solutions we can always 
replace point charges on the ions by rigid charge distri- 
butions without changing any physics if we appropriately 
modify the core interactions as in ea l25l This can be very 
useful because the rapidly varying short ranged parts of 
the Coulomb interaction can often be more accurately 
treated by the same specialized methods used for the 
other strong core interactions, which must be present in 
any realistic model of ionic solutions. 



B. Scaled LMF equation for the OCCHS 

We now apply the general LMF eqElto the OCCHS 
in the special case where the external field (/)(r) = w(r) 
is that resulting from an ion fixed at the origin, given by 
eqO This choice allows us to describe uniform fluids, as 
discussed in Section lllljl We take advantage of spherical 
symmetry and use the Gaussian charge separation of w(r) 
given in easl24l and 1251 

The LMF eq [S] can be naturally rewritten in terms 
of the more slowly varying "perturbation part" of the 
effective field 



'i'Riif) = (\)R{r) ~ uo{r). 



(27) 



If the perfect cancellation argument were exact, then 
4'r{t) — uo{r), or (f>Ri{r) — 0, corresponding to a fixed 
mimic particle at the origin, and the resulting induced 
density in the mimic system would be Apo(''; [uo]) = 
p^/io(r), with ho{r) the pair correlation function in the 
uniform mimic system. The LMF equation corrects this 
approximation by determining a finite short ranged effec- 
tive field perturbation (pjnir), which we can picture as 
arising from a modified solute particle at the origiuf^iii 
that takes better account of the locally averaged effects 
of the slowly varying interactions ui. 

Taking Fourier transforms, and introducing for reasons 
that will soon become apparent a multiplicative param- 
eter a that scales the amplitude of ipRi, a scaled version 
of the LMF equation can be written as 



exp[ 



\{kar]Sn{k), (28) 



where 



SRik) = l + Apo{k; [M) 



(29) 
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For now we simply note that the original unsealed LMF 
equation has a = 1, and view a as a parameter at our dis- 
posal. A similar scaling of the LMF equation for systems 
with short ranged interactions was discussed earlier.^ 

Before giving quantitative results in sections IVl and IVTI 
below, let us discuss some qualitative features of a self- 
consistent solution of eq[2HI Such a solution would deter- 
mine a short ranged effective field, implying a (j)Ri (fc) that 
is finite as fc ^ 0, along with the associated Apo{k; [4>b])- 
By the fundamental assumption of LMF theory in eq |21 
the latter is the LMF approximation to the full A/5(fc; [w] ) 
in the OCCHS. In particular S'/j(&) in eq^iis the LMF 
approximation to <S'(fc) in ea 1161 and when no confusion 
will result, we will simply write S{k). 5ij(fc) should be 
carefully distinguished from So{k) = 1 -f Apo{k; [uq]), 
which equals the structure factor in the uniform mimic 
system with (pm — 0. 



C. Choice of a 

In order that (j)Ri{k) remain finite as fc — > in ea 1281 
it is clear that with any choice of a, the associated Sr{0) 
must vanish identically. (In practice it is not easy to 
ensure this in a self consistent iterative solution of eg 1281 
and we give in the Appendix in ea I A. 31 an alternate but 
equivalent version that is numerically more stable.) Thus 
the LMF theory gives an approximate structure factor 
that always obeys the exact neutrality condition. 
Its expansion at small k has the general form 

Snik) = + B{a)k^ +0{k^), (30) 

resembling eqEI but the coefhcient B{a) of fc^ depends 
on a and does not necessarily obey the exact second mo- 
ment condition. However, by substituting ea I3UI into eq 
1281 we see that the exact result B{a) = k^ is found if a 
is chosen self consistently so that^ 

/3p^<^fl,i(0) = a. (31) 

Thus by proper choice of a, we can guarantee that 
the approximate structure factor S]i{k) also obeys the 
exact second moment condition. We show below in sec- 
tion that with optimal choices of the key parameter 
a, even the unsealed LMF theory with a = 1 often gives 
very good numerical results. However it is conceptually 
important to realize that the LMF approach can be nat- 
urally generalized as in eq |2H1 so that the exact second 
moment condition is satisfied, and this adds essentially 
no numerical costs to the self consistent solution. We use 
ea 1281 along with eq|^in most of the work reported be- 
low, and usually refer to this generalized approach as the 
LMF theory. If we want to emphasize that the second 
moment condition is satisfied, we will refer to the LMF2 
theory, and distinguish this from the original unsealed 
version, which we will call the LMFO theory. 



D. Choice of a 

The ability to choose a a larger than some CTmin allows 
for a consistent application of LMF theory to ionic flu- 
ids, ensuring that the LMF approximation is used only 
for slowly varying parts of the Coulomb interactions. The 
choice of CTmin determines an effective Coulomb core size 
from the core component Wqo (r) of the Coulomb interac- 
tion. This may be larger or smaller than the "physical" 
core size d, which can be varied independently in the 
OCCHS. 

For strong coupling with F ^ 1, we expect consid- 
erable cancellation of the very strong forces from ions 
at distances larger than the nearest neighbor spacing a. 
Thus the effective core size CTmin/a should be of order 
unity and essentially independent of F at large F. 

However, for weak couplings with F <C 1 we would ex- 
pect that any choice of ct > CTmin with CTmin — Fa will be 
sufficient for the LMF equation to give good results. This 
(conservative) choice of CTmin corresponds to the Bjerrum 
length.^ Only on length scales CTmin much less that the av- 
erage neighbor separation a will even the bare Coulomb 
interactions between ions exceed fc^T, which would char- 
acterize an effective Coulomb core size. A more detailed 
argument^^ suggests that we can take even smaller ct, 
including ct = 0, for F <C 1 and still get accurate re- 
sults from the LMF theory. We will see below that these 
qualitative considerations hold true generally. 

In particular, by choosing ct large enough we can guar- 
antee that I3p^4)fa{k) is nonzero only at small wavevec- 
tors, since the Gaussian factor in eq EHl causes rapid 
decay for ^ct > 2. The Gaussian charge distribution 
produces this very efficient localization of [3 (f) iii{k) to 
small wavevectors, and is much superior in this regard to 
most other smooth distributions. 

This property is very important at high density and 
strong coupling where S^^k) can have significant struc- 
ture at ka ~ 27r, where a roughly measures the typical 
distance between nearest neighbor particles. At those 
larger wavevectors that characterize short ranged struc- 
ture in r-space, Pp^ (fiuiik) essentially vanishes for any 
choice of CT > CTmin — «• Thus for such wavevectors we 
have S]i{k) ~ So{k) from a crude Hnear response argu- 
ment. Differences in these functions should show up only 
at small fc, where Sn^k) will satisfy the SL moment con- 
ditions while So{k) remains finite as fc — > 0. However at 
high densities the compressibility in the mimic system is 
small, so that the absolute differences between Sn{k) and 
So{k) remain small even at small wavevectors, as will be 
illustrated in Figure |51 below. 

On taking inverse transforms we then expect that 
ho^r) ~ h{r) holds true to a very good approximation 
at high densities, as was qualitatively suggested by the 
cancellation argument leading to eq El Thus we pre- 
dict a family of uniform mimic systems for different ct > 
CTmin, all of which should give essentially the same short 
ranged structure at high density that closely approxi- 
mates that of the full ionic system. This is a dramatic 
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example showing that the inverse problem of uniquely 
determining the intermolecular potential from /io(?') can 
be ill-conditioned. 

A "molecular-sized" choice of cr ~ Cmin for the mimic 
system is considerably smaller than the typical choices 
made in Ewald sum methods, where a is usually taken 
to be proportional to the simulation system size.^^ Larger 
a values will give equally good results, provided that the 
resulting mimic system is described accurately. However 
generally there is little to gain from such choices, since 
the LMF theory is already consistent for a ~ cr,jii,i, and 
it may be more difficult to treat the longer ranged inter- 
actions in mimic systems with larger a. Thus optimal 
choices for a are generally near CTmin- 

At lower densities and weak coupling there is little 
structure in Sulk) and So{k) at larger wavevectors, and 
we can choose a much smaller (Tmin as argued above 
and still make consistent use of the LMF approximation. 
Very accurate results for SR{k) are again found from a 
self-consistent solution of eq|^ but the long wavelength 
perturbations from /3p^ (j)fii{k) are not damped by low 
compressibility in the mimic system, and ho(r) can dif- 
fer noticeably from h{r). As F — » 0, our theory reduces 
correctly to the exact Debye-Hiickel (DH) limit ii 



OCP: r=6 [LMF2] 




ak 



Figure 2: SR{k) = S{k) at moderate coupling for the OCP 
computed using the MPB theory with different a's and com- 
pared to So{k) for the uniform mimic system. Also shown are 
the associated pf3(j}Ri, which use the scale on the left j/-axis. 
When pl3(j)Ri is taken into account using the MPB theory, 
both choices of a give very similar S{k) that compare well 
with simulation data for the full systemi— 



V. RESULTS AT LOW DENSITY: MIMIC 
POISSON-BOLTZMANN APPROXIMATION 

At low enough bulk densities, the mimic system's re- 
sponse to 4'R{r) can be described using the simple ideal 
gas Boltzmann factor as in eq|Hl so that 

Apo(r;[0fl]) = p^[e-^^«(^)-l]. (32) 

This also represents the LMF prediction for the full sys- 
tem's p^h[r), and requires only that second and higher 
order virial corrections to the mimic system's pair cor- 
relation function can be ignored. When eq is sub- 
stituted in the LMF equation 1281 a closed equation for 
4>R{r) results. A self consistent solution is readily found 
by iteration, using the equivalent but numerically more 
stable version of the LMF equation in eg IA.3I 

Eauation l32l is the same structural approximation that 
is used in the Poisson-Boltzmann (PB) theory.^ Indeed if 
we set a = 1 and cr = in the LMF eauationl28lcombined 
with eg 1321 the results reduce exactly to those of the stan- 
dard PB theory. The PB theory thus results from taking 
the full Coulomb interaction of eq as the perturba- 
tion ui in the LMF equation and using the Boltzmann 
approximation for the density response. 

We refer to the low density limit of our theory, where 
eg 1321 is used in eg 1281 as the mimic Poisson-Boltzmann 
(MPB) theory. The MPB theory differs from the PB the- 
ory only by the choice of a yielding a consistent mimic 
system along with a choice of a that ensures that the sec- 
ond moment condition is exactly satisfied. As we will see, 
these simple modifications greatly improve the accuracy 
and range of validity of the MPB theory. 



A. MPB theory for OCP 

Consider first the OCP, where there is no length scale 
in the potential to suggest an intrinsic core size. A qual- 
itative discussion of the choice of the effective size CTmin 
was given above. In practice it is easy to determine CTmin 
by solving eg 1281 using successively larger values of a. For 
CT < CTmin the rcsults vary strongly with ct and are gen- 
erally very inaccurate. But for all ct > CTmin the LMF 
theory is consistent and should give very similar predic- 
tions for the full system's structure as exhibited in S'fl(fc), 
even though the effective fields (}p^(j)Ri{k) and the uni- 
form mimic systems' structure factors 5*0 (fc) can still vary 
strongly with ct. This is illustrated in Figure |21 for the 
state with F = 6, where the convergence of the results 
for a/a — 0.69 and a /a ~ 0.73 is shown. 

The LMF equation itself would continue to give (even 
more) accurate results for larger ct if an accurate the- 
ory for the structure Apo{r; [0^]) induced by a given (f>ii 
is used. However the simple Boltzmann approximation 
used in eo l32l for this quantity must fail at higher densities 
where there are significant correlations between mimic 
system particles. This sets a CTmax above which the re- 
sults of the MPB theory become inaccurate, very roughly 
estimated by P^ctJ^j^^ ^0.1 as for hard sphere fluids. 

As F increases, eventually the CTmin needed for the ac- 
curacy of the LMF approximation exceeds this CTmax and 
the MPB theory fails. The internal consistency or in- 
consistency as CT is varied is very evident from the MPB 
theory itself. In practice we find very good and consis- 
tent results for all F < 6 and the slight differences in 
the S{k) curves in Figure |21 and deteriorating results at 
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0CP:r=6 [o/a=0.69] 




Figure 3: Moderate coupling OCP structure computed us- 
ing the MPB theory. In the left graph, the LMF2 S{k) is 
compared with the result of the generalized Debye Hiickel 
theoryj-^i (GDH) and the usual PB theory. The right graph 
makes the same comparison for g{r) and also shows the effec- 
tive field perturbation fScfimir), which uses the scale on the 
left j/-axis. Both the LMF2 and GDH solutions satisfy the 
second moment condition while PB does not. The GDH re- 
sult is expressed as an expansion and computed up to its Z = 6 
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Figure 4: Low density OCCHS structure. The LMF2 g{r) 
from the MPB theory is compared with the result of density 
functional theory^^ (DFT) and the PB and DH approxima- 
tions. Though constrained to be zero inside the hard core 
by a boundary condition, the DH g{r) has a negative region 
near r = d, where uo(r) is also rapidly rising, with a value 
of —1.62 at contact. The PB and DH theories fail to capture 
qualitatively the onset of oscillation in g{r) at this moderate 
coupling strength. 



larger a indicate that we are near the upper limit of T 
where the MPB theory can be trusted. This represents a 
surprisingly strong coupling, since the lowest order Boltz- 
mann approximation for the structure in eq 1321 is used, 
and shows the virtues of choosing a mimic system. 

As shown in Figure 13 the results of the MPB the- 
ory for r = 6 are in very good agreement with com- 
puter simulations^^ and compare very favorably to those 
of the usual (nonlinear) PB theory or the generalized 
Debye-Hiickel (GDH) theory developed by Levin and 
coworkersi^ Note that the MPB theory, unlike the usual 
PB approximation, can predict oscillations in both S{k) 
and g{r) from the self consistent determination of 4iRi{r) 
despite using only the lowest structural approximation, 
eqOa 



B. MPB theory for OCCHS 

We now turn to the OCCHS. This has more complex 
structure because of possible competition between corre- 
lations induced by the hard cores and the soft repulsive 
Coulomb interactions. We follow the usual convention 
where lengths are measured in terms of the hard core 
diameter d. 

In Figure 0] we compare the results of the MPB the- 
ory to MC simulations, to results of a density functional 
treatment, and to the PB and DH theories for a low 
density state with ry = 0.02 but with a moderate ionic 
strength F = 5.43. The MPB theory gives excellent re- 



sults with a molecular-sized a = l.Sd, noticeably better 
than those of the considerably more complicated density 
functional theory,^^ 

The right panel in Figure [S] illustrates the role of a for 
this state, and shows that the self consistent choice of 
a = 1.2 can ensure that the exact second moment con- 
dition is satisfied, though on the scale of the graph the 
differences between the Sji{k) with a = 1 are hardly vis- 
ible. The left panels show that for much weaker coupling 
with F = 0.54 even the usual PB approximation (with 
cr = 0) gives good results, indicating that at weak cou- 
pling the choice of cr is not so important. However, unlike 
the MPB theory, the PB approximation can satisfy the 
second moment condition exactly only in the limit F ^ 0, 
where it reduces to the DH approximationji 

VI. RESULTS AT HIGH DENSITY AND 
STRONG COUPLING 

A. Structure in Mimic and Full Systems 

At high densities we expect that the local structure 
in r-space of even the uniform mimic system, where 
all corrections from are ignored, will closely resem- 
ble that of the full system as suggested in eq This 
is illustrated in Figure IHl for the OCCHS for high den- 
sity states with moderately strong couplings. Canonical 
Monte Carlo simulations are carried out to obtain the 
uniform mimic system's correlation function go(r)^ and 
these are compared to previous simulations for the full 
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Figure 5: The upper and lower left graphs show the MPB 
and PB approximations for the OCCHS g{r) compared to 
MC data-- for the full system for weak and moderate ionic 
strengths. The PB approximation is satisfactory only at weak 
couplings. The right graph shows that varying a so that 
the second moment condition is satisfied in the MPB the- 
ory changes the amplitude of (5p(j)Ri{0) but in this case the 
effects on S{k) are hardly visible on the scale of the graph. 



systemjS where Ewald sum methods were used to ac- 
count for the long ranged interactions. Because of the 
short ranged interactions uq in the mimic system, our 
simulations are completely straightforward and no Ewald 
sums or other special treatment of the periodic boundary 
conditions are required. Again a molecular-sized choice 
of a = d of order the nearest neighbor spacing suffices. 

Also shown on the same graphs are the bare ion poten- 
tial (3w{r) and the mimic potential (3uQ{r). Despite the 
much smaller amplitude of the latter and its much shorter 
range, the mimic go{r) has a striking resemblance to the 
full system's g{r), and the differences are hardly visible 
on the scale of the graph. 

At the strongest couplings, the correlation functions at 
both densities have a first peak shifted away from con- 
tact with the embedded hard sphere. Such a correlation 
function is very different from the correlation function 
of a hard sphere fluid, which has its maximum at con- 
tact, and shows that the strong short ranged parts of 
the Coulomb interactions can compete with packing ef- 
fects from hard cores even at high density. This also 
emphasizes the importance of having the softer piece 
Wqa(r) = q^ericir / a) / er outside the hard core in our 
mimic system potential (3uq (r) in eq 1251 in order to re- 
produce correlation functions in the OCCHS, especially 
for strong coupling states. 

There have been several previous empirical attempts 
to fit correlation functions for Coulomb systems at high 
density using effective short ranged systems. The DH 
limit might suggest that a generalized Yukawa fluid could 
be usefulfS^ but the results were not very accurate and 
there was no systematic way to determine parameters for 
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Figure 6: OCCHS correlation functions g{r) at stronger cou- 
pling strengths T and large packing fractions for the full 
and mimic systems as determined by MC simulations. Note 
the maximum in g{r) away from contact for F — 70, indicat- 
ing the strength of the Coulomb repulsions. Also shown are 
I3w{r), the full potential of the ionic fluid, and the mimic po- 
tential l3uo (r) , both of which use the scale on the right j/-axis. 
There is a hard core interaction for r < d. 



the effective potentials. 

Most relevant to our work are ion reaction field (RF) 
methods, where an effective finite-ranged interaction 
w^Q^ {r) was originally determined from the electrostatic 
potential of a positive point charge surrounded by a neu- 
tralizing uniform spherical charge distribution with a ra- 
dius rc.-22. Good results for correlation functions were 
found using the RF method in several applications at 
high density, though some spurious oscillations were seen 
near the cutoff rc- These were attributed to the dis- 
continuity of the second and higher order derivatives 
of w^Q^ {r) at Tc and better results were found using a 
smoother "charge cloud" distribution that had disconti- 
nuities only in fourth and higher order derivatives at the 

cutoff, 

Our WqQ in eg 1261 can be similarly interpreted in terms 
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OCCHS: r=70, Ti=0.4 
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Figure 7: Illustration of mimic system behavior. OCCHS 
correlation function for the state F = 70 and 77 = 0.4 com- 
pared to those of two different mimic systems with different 
a values. Also shown are (3w{r), the full potential of the ionic 
fluid, and the mimic potentials Puo{r), all of which use the 
scale on the right y-axis. 

of the potential arising from a positive point charge 
surrounded by a neutralizing Gaussian charge distribu- 
tion. All derivatives of Wqo are continuous because of the 
smooth cutoff, and by construction the associated per- 
turbation Wqi decays very rapidly in fc-space. It is the 
latter property that fundamentally leads to mimic system 
behavior with a proper choice of a. Our work thus pro- 
vides a conceptual framework for understanding why RF 
methods can work as well as they do in some cases and 
how results can be significantly improved, especially at 
lower densities or in nonuniform environments by using 
the LMF theory. 

Figure \7\ gives a more detailed comparison of the 
structure of the high density/strong coupling state with 
r = 70 and 77 = 0.4 to that of two different mimic sys- 
tems with a = d and a = l.bd. Despite the fact that 
the (repulsive) potential of the latter is always greater 
than or equal to that of the former, both mimic systems 
have very similar correlation functions that agree very 
well with that of the full system, which can be viewed 
as the limit ct = 00. Thus for high density states all 
mimic systems with a > CTmin have essentially the same 
short-ranged structure in r-space. 

Figure |S1 gives similar results for the OCP at very 
strong coupling strengths F = 80 and F = 140. We find 
excellent agreement with simulations of the full system^i 
using a mimic system with a/a = 1.4. As we would ex- 
pect, for small enough a the good agreement fails, as 
illustrated by results for a /a — 1. 
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Figure 8: High coupling strength OCP correlations for the 
full and mimic systems. For a = 1.4a, go{r) is essentially in- 
distinguishable from the full system's g(r). However a smaller 
a = a fails to mimic the full system's correlations. Note that 
around F ~ 170 the OCP starts to freeze.^ Our simulations 
indicate that the mimic system with a — 1.4a also freezes at 
around the same F. 

Figure compares the mimic structure factor So{k) 
and a simple estimate for SR{k) — S{k) based on a linear 
response treatment of the effects of (3p^ 4>Bi{k)- Only 
at very small k as revealed in the inset can any differences 
be seen. The linear response treatment turns out give an 
Sn.{k) that satisfies exactly both the zeroth and second 
moment conditions with a = 1, and the results converge 
to So{k) at larger k controlled by the factor exp[—j{ka)'^] 
arising from our choice of a mimic system. These fea- 
tures would be found in any more exact treatment and 
suffice at high densities to give a very accurate estimate 
of SR{k). This also suggests that the results from other 
approximate theories may be improved by use of a good 
mimic system. Indeed, as already shown in Figure El 
the simplest possible theory where the effects from (jjjn 
are ignored completely, already gives excellent results for 
short ranged correlations in r-space. 



B. Internal Energy at High Density 

With accurate approximations for S{k) and h(r) in 
hand, it is straightforward to calculate thermodynamic 
properties by integration. The simplest of these is the 
internal energy, given by eq 1191 One can always solve 
the LMF equation to obtain h{r) and use this integral to 
compute the internal energy of an ionic system, and we 
would expect very accurate results. Here we show that 
because of the great similarity of the local structure of 
the mimic and full systems in r- and /c-space at high den- 
sity, we can obtain an accurate estimate of the energy in 
terms of the mimic system's energy and a simple analytic 
correction without explicitly solving the LMF equation. 
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OCCHS: r=20, ri=0.343 [LMF+linear repsonse] 




Figure 9: Linear response theory is used to approximate the 
change in SR{k) = S{k) in the mimic system induced by the 
the perturbation 0iji(fc). With a — 1, the linear response 
theory for S{k) satisfies both SL moment conditions. The 
inset is a blown-up view of the structure factor at small k 
where the differences in So{k) and S{k) can be seen. 



Equation ^1 can be rewritten in fc-space and the 
Coulomb interaction separated into mimic and pertur- 
bation parts, so that 



1 1 



N 



2 (27r)3 
1 1 



2 (27r)s 



-i^'^^^'p^'hik). (33) 



In the first term of eq|22l because p^ho{k) differs from 
p^h{k) only at small k, where the factor 1 — e^i^'"^^ also 
approaches zero, we can replace the latter by the former 
with little error. In the second term, e"*^'''^-' decays 
very quickly at higher fc, so only the small-fc features of 
p^h{k) ~ — l+fc^/A;|, from the SL moment conditions are 
needed to have an accurate estimate of the integration. 

The first term then gives the internal energy of the 
mimic system 



(3EI 



ex 



N 



drua{r)p^ho{r), 



(34) 



(with a background contribution —Tr/3q'^a'^p^/e arising 
from the use of Hq in the integral rather than go), while 
the second term corrects the mimic system's energy and 
can be integrated analytically to give 



-1 



N e^cr V ^ ' (fcD^T)^ 
Thus the internal energy can be estimated as 

PE^'' PEf"" 



N 



N 



N 



(35) 



(36) 



r 




a/a 


a/d Full MC Theory [eq 




20 


- 


1.4 


- 


-16.67 


-16.59 


-8.66 


20 


0.343 


1.4 


1 


-17.17 


-17.12 


-9.20 


20 


0.4 


1.4736 


1 


-17.33 


-17.27 


-9.73 


70 


- 


1.4 


- 


-60.81 


-60.72 


-32.65 


70 


0.343 


1.4 


1 


-61.09 


-61.02 


-32.94 


70 


0.4 


1.4736 


1 


-61.32 


-61.24 


-34.56 


80 




1.4 




-69.69 


-69.64 


-37.54 


125 




1.4 




-109.73 


-109.74 


-59.50 


140 




1.4 




-123.09 


-123.13 


-66.85 


160 




1.4 




-141.72* 


-141.57* 


-77.23 



Table I: Excess internal energy BE^'' /N. The full MC data 
are taken from referencesl^^ andl3ll At F = 160, both the full 
and the mimic systems are near solidification, and the results 
depend on initial conditions. 



Results for this approximation are compared to MC re- 
sults in Table U for a variety of high density/strong cou- 
pling states. The Gaussian charge distribution is the key 
to the accuracy of eql^HI both because of its fast decay 
in /c-space and its use in revealing an excellent mimic 
system for the full ionic system. 

Note that at high densities or ionic strengths, the De- 
bye wave vector ku in eq El can be very large, mak- 
ing neutrality the major contribution to PEf^/N ~ 
—(3q^/ {e^/ncj) . Though the physical interpretation of this 
term is very different, its limiting value is the same as the 
self-interaction correction in the Ewald sum methodic A 
similar correction was used in ion RF methods 

At high density where the compressibility xt of the 
mimic system is very small, we can obtain a good esti- 
mate of the internal energy by simply replacing h{r) by 
/io(r) in eq^l where because of the background subtrac- 
tion, a finite result is found for any short ranged ho{r). 
Separating the Coulomb interaction as in eg 1331 and us- 
ing only the k = result p^ho{0) = — 1 + ksTpxT in the 
second integral we find analogous to eg 1351 



N 



(-1 + ksTpxT) 



(37) 



which agrees with eg 1351 in the relevant limit of large k^ 
and small xt- 



VII. FINAL REMARKS 

It is straightforward to apply these idea to charge and 
size asymmetric primitive models^^ or to "simple molten 
salt" models^ with softer repulsive cores. We have also 
used the LMF approximation to look at the OCP near a 
charged hard walli^ Here very accurate results are found 
both in the weak and strong coupling limits, and the 
more complicated pair level theory recently introduced^^ 
is not required. Details will be presented elsewhere, and 
we only make a few general remarks here. 
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We have found for size asymmetric primitive models 
that the simplest choice of a single tr parameter for all 
species gives excellent results. This can be understood 
as a consequence of Stillinger and Lovett's fundamen- 
tal insight that for general ionic mixtures universal con- 
sequences of the long ranged Coulomb interactions can 
be seen in the small wavevector behavior of the charge- 
charge correlation function. By using the smeared charge 
distributions implied by a proper choice of u, we arrive 
at a smeared charge-charge correlation function that has 
significant structure only at small wave vectors. The 
LMF theory can then reproduce the exact long wave- 
length behavior found by SL, and the slowly varying 
smeared Coulomb perturbations have little effect on the 
shorter wavelength correlations induced by the modified 
ion cores in the mimic system, just as illustrated here for 
the OCCHS. 

The main complication that arises for such mixtures is 
that the resulting mimic systems will have short ranged 
attractive as well as repulsive interactions, but this is 
required if we want the mimic system's structure to re- 
semble that of the full Coulomb system. However, the 
LMF theory can then to a very good approximation take 
care of the "universal" long ranged parts of the Coulomb 
interactions, which cause major conceptual and compu- 
tational problems in most approaches, while leaving a 
nonuniversal, but surprisingly short ranged problem to 
be treated by whatever means are available. Simulations 
are straightforward and some recent theoretical develop- 
ments for treating systems with strong but short ranged 
interactions look very promisingi2& We believe the ideas 
presented here offer a powerful general perspective, and 
are actively pursuing their consequences for static and 
dynamic properties of fluids with both short and long 
ranged interactions. 
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Appendix: NUMERICALLY STABLE VERSION 
OF LMF EQUATION 

The LMF equation 1281 will produce the desired short 
ranged 4'Ri{r) with a finite value of (j}Ri{Q) only if Sii{k) 
rigorously vanishes at fc = as in eql^OI which of course 
is the exact result. However such self consistent equa- 
tions are usually solved by iteration and any small errors 
in an intermediate approximation to Sfi{k) at small k are 
greatly amplified. This can lead to numerical instabili- 
ties. The following rewriting of eq|2H|can alleviate this 
problem. We can remove the sensitivity at small k by 
multiplying both sides of eq|2H|by k'^, giving 

k''l3p''$Ri{k) = aklexp[-l{kaf]SR{k). (A.l) 

This equation suffices to determine (fc) everywhere 
except near fc = 0, where (f>Ri is assumed to be regular. 
We next formally write an identity involving (jjRi (k) that 
remains finite as /c — + 0: 

K^Pp''$Riik)^K^Pp''$Ri{k), (A.2) 

where if is a (real) constant wavevector. (More generally, 
we can multiply both sides by a known real function of 
k that does not vanish as /c ^ 0.) We now add these two 
equations as written and divide by k^ + K'^, which yields 
our final result 

f3p^Mk) = -^^eM-\{k^f]SR{k) 

+ ^^^P''Mk). (A.3) 

This equation has no divergences at small k and is stable 
when iterated, with a choice of K of order ki). A con- 
verged solution will produce a S'/?,(fc) that vanishes iden- 
tically at fc = with Pp^ipRiiO) finite. From eg 1311 this 
quantity is in fact of order unity when the second moment 
condition is satisfied. This way of rewriting equations of 
this kind was suggested to us by Kirill Katsov. 
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